The Bethe ansatz, which was introduced in 1931 by Hans Bethe, has become a powerful method to obtain exact solutions of one-dimensional (1D) quantum many-body systems. In 1963, Lieb and Liniger^[1] solved the 1D many-particle problem of δ-function interacting bosons by the Bethe’s hypothesis. The ground state, the momentum, and the elementary excitations were obtained for this model by using the Lieb–Liniger solution. In this context, a significant step was made on the discovery of the grand canonical description of this Lieb–Liniger model by Yang and Yang in 1969.^[2] Now, this grand canonical approach is called Yang–Yang thermodynamic method. The Yang–Yang thermodynamics of the Lieb–Liniger Bose gas provides benchmark understanding of quantum statistics, thermodynamics, and quantum critical phenomena in many-body physics, see a review.^[3,4] In the context of ultracold atoms, the 1D Bose gas with a repulsive short-range interaction characterized by a tunable coupling constant exhibits rich many-body properties. This model thus becomes an ideal test ground to explore fundamental many-body phenomena ranging from equilibrium to nonequilibrium physics in the experiment.^[5–10]

Despite Lieb–Liniger is arguably the simplest integrable model, the calculation of the correlation functions is extremely challenging due to the complexity of the Bethe ansatz many-body wave function of the model. The study of correlation functions has long been being an important theme in the physics of ultracold quantum gases since they provide information of quantum many-particle interference and coherence beyond the solely spectra of the systems. Therefore, there has been much theoretical and experimental interest in the local, non-local, and dynamical correlation functions at zero and finite temperatures via numerous methods based on exactly solvable models, see Refs. [11]–[22].

For sufficiently strong interaction and sufficiently low density, the 1D Lieb–Liniger gas enters the Tonks–Girardeau (TG) regime, in which bosons behave like impenetrable particles (hard-core bosons). Such impenetrable bosons behave mostly like the free fermions that build up the Girardeau’s Bose–Fermi mapping.^[23] In fact, the 1D Lieb–Liniger Bose gas with the interacting strength c_B can map onto the fully-polarized fermions with a p-wave interaction of strength c_F = 1/c_B.^[24] As a result of the Bose–Fermi mapping, the energy spectra of the Bose and corresponding Fermi systems are identical at the TG regime. The observables that can be given in terms of the local density are identical for both systems, such as dynamical density–density correlation function, see Ref. [25]. However, this mapping for the off-diagonal correlation function does not like to be true, for example, the momentum distribution. The Bose–Fermi mapping has tremendous applications in the study of strongly interacting quantum gases of ultracold atoms. An expansion of the reduced density matrices for 1D bosons has been first derived by Lenard^[26,27] using Bose–Fermi mapping at Tonks Girardeau regime. After that, there were superb findings in the literature.^[28–35]

In the context of correlation functions, the Fourier coefficients have an important physical significance. By taking Fourier transform of the one-and two-particle correlation functions, the momentum distribution and the static structure factor of the system can be obtained, respectively. The first few terms of the short-range series expansion of the one-particle correlation function have been calculated in Refs. [28], [36], and [37], where c₁ = 0, c₂ = 1/2(e(γ) − γe′(γ)), c₃ = γ²e^′(γ)/12, etc. A further study on the connection between the non-local one-body and local three-body correlation functions was presented in Ref. [38]. The systems described by the Hamiltonian with short distance interaction (δ-function type) show singular behavior of the wave function when two or more colliding particles coincide. Consequently, the large momentum tail of the momentum distribution has a universal power-law decay, i.e., w(p)∼C/p⁴, where C is an extensive quantity called Tan’s contact.^[39–41] Beside the power law tail of the momentum distribution, the contact characterizes many universal properties of the systems that are independent of the details of the contact interaction and the dimensionality,^[42] such as the energy adiabatic relation, the universal connection to the structure factor, and the pressure relation. The contact measuring the correlation of such systems of this kind in short distance limit is proportional to the probability of finding two atoms at zero separation,^[39–41] see more studies on the contact in Refs. [43]–[45]. The contact can be computed from the free energy of the system via the Helmann–Feynmann theorem,^[46] or from the prefactor of the high-frequency power-law decay of the single particle correlator,^[47] or by the operator product expansion field theory method.^[48] Very recently, p-wave contacts were found to exhibit universal relations in the ultracold atomic systems with a p-wave interaction.^[49–51] In this scenario, the momentum distribution, s-wave and p-wave contact relations have been attracted a great deal of attention from theory and experiment.^{[36,37,42,46,52–56]}

On the other hand, the static structure factor contains information about the static correlation properties of a system and directly relates to the pair distribution function or the normalized density–density correlator or spin–spin correlator.^[25] In the low-momentum region and the infinite repulsive regime, hydrodynamic theory predicts that the static structure factor should behave linearly with momentum k.^[42] The static structure factor at small momentum k is related to the sound velocity (v_s) by S(k) = ħ|k|/(2mv_s).^{[12,21,22,25,42,57]} For high momentum, the structure factor always converges to 1.^[42] The Fourier transform of the time-dependent density–density correlation function is known as the dynamical structure factor,^[21,58,59] and expresses the probability of excitation by an infinitesimal perturbation.

The duality between interacting spinless bosons and fermions in 1D^[24,25,32] can be adopted to the framework of the generalized fractional statistics, which interpolate continuously between bosons and fermions. In fact, while three dimensions particles can be only bosons or fermions, in lower dimensionality it has been recognized that particles with fractional statistics intermediate between bosons and fermions^[60–66] due to the transmission between dynamical interaction and statistical interaction. The 1D strongly repulsive Bose gas at low temperatures is equivalent to a gas of ideal particles obeying the non-mutual generalized exclusion statistics (GES) with a statistical parameter α = 1 − 2/γ, where γ is the dimensionless interaction strength.^[66] Using such a mapping, the quasimomenta of N strongly interacting bosons map to the momenta of N free fermions via with i = 1,…,N.^[66] By using the fractional statistical parameter α, the M-body local correlation functions and non-local correlation functions in terms of the wave function of M bosons at zero collision energy and zero total momentum have been computed analytically up to next-to-leading order.^[19] In this paper, we analytically calculate various correlation functions in terms of the statistical parameter α. In particular, we calculate the two-point correlation function, the power-law tail of the momentum distribution, and discuss the structure factor of the 1D Lie–Liniger model of ultracold atoms with finite strong interaction strength. For Tonks gas, i.e., α = 1, our results are equivalent to those of the Tonks–Girardeau gas.

The paper is organized as follows. In Sections 1 and 2, we briefly introduce the Lieb–Liniger model and recall some results of higher-order local correlations for the 1D interacting bosons. In Section 3, we present the derivation of the nonlocal 2M-point correlation function. In Section 4, we prove the power-law decay of the momentum distribution of the 1D interacting bosons at large momenta. The large momentum tail of the momentum distribution determines Tan’s contact in term of the fractional statistical parameter α. In this section, we also discuss the static structure factor. The last section presents our conclusion.

A model of 1D spinless N bosons with δ-function interaction is described by the Hamiltonian

With the periodic boundary condition ψ(0,x₂,…,x_N) = ψ(x₂,…,x_N,L), in the domain 0 ≤ x₁ ≤ x₂ ≤ ··· ≤ x_N ≤ L, the Betha ansatz wave function can be written as^[1]

For a finite large interaction strength, the quasimomenta of the bosons k₁, k₂,…,k_N deviate from pure Fermi statistics, they obey the non-mutual general exclusive statistics (GES).^[60,61,66] The deviation from Fermi statistics is determined by the non-mutual GES parameter α = 1−2/γ.^[66] If the total momentum k₁ +··· + k_N = 0, we have the following relation:

We first would like to recall some results of local M-body correlation functions which have been studied in Refs. [11]–[22]. In particular, we discuss the results of the correlation functions which were recently derived by using the GES relation Eq. (4). It turns out that the GES is convenient to calculate the local M-body correlation functions.

3.1. Strong coupling regime

The local M-particle correlation is a measure of the probability of observing M particles at the same position and defined by

In the strong coupling limit, one can expand the BA wave function Eq. (3) in terms of the power of 1/c, from which the numerator and denominator can be calculated analytically^[19] Here, is the Vandermonde determinant and , for all x’s. In this case, using the GES Eq. (4), , and making a scaling change for i = 1,…,N, thus we can rewrite Eq. (7) as to the order O(c^{−M(M − 1)−2)}. Here, . Since is a Slater determinant, by applying Wick’s theorem and changing variable k^F = 2π n z, then the M-body correlation function reads where the dimensionless parameter γ = c/n, the liner density is given by n = N/L and , T and μ are the temperature and the chemical potential of the 1D Bose gas, respectively. Calculating the multiple integrals in Eq. (9) by using the orthogonal polynomials in the random matrix theory,^[67] we can give the M-particle local correlation function as where h_i’s are the norm-squares of the monic orthogonal polynomials satisfying . Considering the distribution of particles at the ground state and applying the Sommerfeld expansion, we have general formulas for the higher-order local correlation function in two temperature regimes^[19] where T_d ≡ ħ²n²/(2mk_B) is the quantum degeneracy temperature. In the above calculation, the statistical parameter α in the quasi-momenta of the hardcore bosons can be extracted as a factor giving the next-order corrections in the correlation function (11). In such a way, the integrant in the correlation function (8) reduces to the Vandermonde determinant with respect to the the wave function of the free fermions. Then the Wick’s theorem and random matrix theory are used to obtain the explicit form of M-body correlation function (11). However, it would be extremely difficult to obtain higher-order correction terms beyond the form of Eq. (11).

3.2. Weak coupling regime

For the weak coupling case, one can easily use the wave function (3) to calculate the M-body correlation function

So, the M-particle local-correlation function becomes where .

In fact, the Bogoliubov method is very convenient and useful for the study of the weak coupling Bose gas in any dimension. In this Bogoliubov approach, one considers a large fraction of particles in the lowest quantum state and excites a small fraction of particles out of the Bose-condensed particles with weakly repulsive interacting Bose gas. Therefore, in this regime, the field operator can be considered as a sum of the condensate part ψ₀ and a non-condensate part ψ′, such as . The density of the system is n = n₀ + n′, where n₀ = ⟨ψ₀ ψ₀^*⟩ and n′ = ⟨ψ′^†ψ′⟩ represent the condensed particle density and non-condensed particle density, respectively. For a uniform condensate, the ground state solution is , where the chemical potential is given by the relation μ = n₀g, see a detailed study given in Ref. [13]. In the weak coupling regime, γ ≪ 1, the mean field interaction energy per particle is proportional to ng, where , and it is proper to define the correlation length . At temperature T≪ T_d, over a wide range of parameters l_c ≪ l_ϕ, we can have the condition .^[13] Here, l_ϕ is the phase coherence length and T_d = ħ²n²/2m is the quantum degeneracy temperature. In this case, the condensate fraction contains the contribution of excitations with momenta k ≤ k₀ ≪ l_c⁻¹ and non-condensate parts provided excitations with k ∼ l_c.^[13,15] In order to complete our discussion on the local M-body correlation, the necessary calculation of the correlation function is shown in Appendix B. What below is a summary of the correlation functions based on the Bogoliubov approach.

The M-particle correlation function for 1D bosons at weak coupling regime is given by

dealing with the M^th order density n^M = n^M − Mn^{M − 1} ⟨ψ′†ψ′⟩ + O(⟨ψ′†ψ′⟩)². Equation (B5) shows that depends on the normal density ⟨ψ′†ψ′⟩ and anomalous density ⟨ψ′ψ′⟩ of the non-condensed part. Considering the form of the non-condensate part, , can be converted into by Bogoliubov transformation, where are operators of excitations which obey the usual commutation relations, ε_k are their eigen energies, and u_k, v_k are their eigenfunctions. This result illustrates that depends on the occupation number of the excitation N_k, excitation energy ε_k, and system energy E_k. In fact, varies with the temperature range since N_k = 0, at zero temperature and N_k ≠ 0, ε_k = ħkv_s at non-zero temperature. Here, is the sound velocity and the distribution N_k takes the form (e^ε_k/T − 1)⁻¹ for the temperature T ≪ μ and T ≫ μ, respectively. However, E_k = ħ²k²/2m and μ ≈ n₀g at low temperature.

This study follows the method developed by Gangardt and Shlyapnikov in Ref. [13]. The elementary excitation of the Bose-condensed system involves vacuum fluctuations and thermal fluctuations. At zero temperature, thermal fluctuations are suppressed and vacuum fluctuations dominate. This scenario becomes reversed at . However, both vacuum and thermal fluctuations create excitations at T ≪ μ.

According to these configurations, and after a length algebra, the M-particle correlation, Eq. (15), is given by

These results are not valid for , belonging to quantum decoherence regime. The detailed calculation is presented in Appendix B. This result has been reported in Ref. [13]. The study of 2-body correlations was also reported in Ref. [14], while the 2-body, 3-body correlation functions at zero temperature were studied in Ref. [15].

The M-particle non-local correlation is a measure of the probability of observing M particles at different points and defined by

4.1. Strong coupling regime

We first consider the non-local correlation function of strongly interacting 1D bosons with c = ∞. Let us shorten our notation g_M(x₁,…,x_M) = g_M(x₁,…,x_M; 0,…,0). According to the Girardeau’s Bose–Fermi mapping,^[23] we have in the limit c → ∞. Thus the M-particle correlation is rewritten as

where N-particle fermion wave function has a determinant form and the norm of the wave function can be expressed as and the normalization factor . By changing the variable as , equation (18) becomes

We further observe that the elements of the determinant in Eq. (19) satisfy

and A = N is a constant. Therefore, according to the theorem, the M-partical non-local correlation function can be written as which can be used for calculating different nonlocal correlation functions of Tonks–Girardeau Bose gas.

4.2. Subleading contribution in 2-body correlation g₂

Forrester et al.^[30] have derived the O(1/c) order correction to the two-body non-local correlation function by expanding the momentum of bosons in the 1D Lieb–Liniger model with a finite strong repulsion. This subleading contribution can also be retrieved by using the fractional statistics following their analytical process. When we consider the non-local correlations, the positions of the particles (i.e., the permutation orders in the wave function) are very important. In the local correlation function, we have used the GES, , the coordinates of the particles are demanded by a rescaling . In contrast, such a rescaling of the particles’ coordinates in a calculation of the non-local correlation function works only for the strong repulsion limit, i.e., x_i/c ≪ 1. For calculating the non-local correlation function, we will use Forrester et al.^[30] expansion method in term of the statistical parameter α.

For our convenience, we fix the number of total particles N + 2 in the domain 0 ≤ y ≤ x₁ ≤ … ≤ x_j ≤ x ≤ x_{j + 1} ≤ ··· ≤ x_N ≤ L and label the particle positions as . Then the two-body non-local correlation function reads

For our convenience in notation, we can denote the wave function in Eq. (5) as Thus, the operators in the wave function can be expanded as

Let y = 0, and depends only on x₁,…,x_N variables. At the ground state, the system has zero total momentum, i.e., . The wave function of free fermions vanishes when two particles coincide, i.e., at x_l = 0,L,x_j(j ≠ l). Therefore,

Taking together these conditions with Eq. (29), we may show that and where .

Furthermore, the coefficient of in Eq. (32) can be simplified in terms of the GES parameter α. Equation (32) has been derived in terms of 1/c corrections.^[30] A detailed calculation of the integral term (in Appendix C) gives

where and are the two-body and three-body non-local correlation functions of N + 1 free fermions, which have been given in Eq. (25). Finally, equation (32) can be viewed as the two-body non-local correlation function of N particles, namely, We observe that the two-body correlation function essentially depends on the two-body and higher order non-local correlation functions of the free fermions and their derivatives, which provide insights into the many-body correlations in strongly interacting bosons.

Various methods can be used to calculate the contact,^{[39–41,43–45]} for example, it can be obtained from the free energy or ground state energy via the Helmann–Feynmann theorem,^[46] or from the single particle correlation function,^[47] from the operator product expansion field theory method,^[48] etc.

5.1. Large momentum tail of the momentum distribution

The single-particle momentum distribution w(p) at momentum p = ħ k is the Fourier transform of the one-particle density matrix ρ(x_j,x_k). In order to build up a connection between the Tan’s contact and the power law decay with distance in the large momentum tail of the momentum distribution, we first calculate the wave function in momentum space. To this end, we introduce the center of mass coordinate R = (x_j + x_k)/2 and relative coordinate x_kj = x_k − x_j of the j,k pair of particles with other fixed N − 2 coordinates {x_{ii ≠ j,k}. The k^th particle in the neighbor of the j^th particle, such as x₁ ≤ x₂ ≤ x₃ ≤ … ≤ x_{k − 1} ≤ x_j ≤ x_k ≤ x_{k + 1} ≤…≤ x_N and x₁ ≤ x₂ ≤ x₃ ≤…≤ x_{k − 1} ≤ x_k ≤ x_j ≤ x_{k + 1} ≤…≤ x_N, the wave function for 1D strongly interacting bosons in Eq. (3) can be approximately written as

Following symmetric conditions, we can present Eq. (35) in general forms as

for bosons in the 1D Lieb–Liniger model, where sgn(x) is +1 for x > 0 and −1 for x < 0, ψ that corresponds to the wave function Eq. (3) at .

The Fourier transform of the one-particle correlation function is the momentum distribution

where p = (2π/L)s and s is an integer. This is equivalent to calculating first the many-particle wave function in momentum space for one atom and then integrating out the rest of the coordinates from its square modulus. In order to calculate the asymptotic behavior (the large p tail) of the momentum distribution,^[42] we first consider the normalized wave function with respect to the first particle in momentum space where p = (2πħ/L)l with l an integer. For a repulsive interaction, the diagonal terms become much larger than the off-diagonal terms due to the periodic boundary condition and symmetry. For a periodic function F₁(x₁,…,x_N);(x_k − x₁) sgn(x_k − x₁), defined on the interval [0,L], where is a regular function, we thus rewrite the square modulus as After a length calculation, the asymptotic behavior of the momentum distribution is given by

Applying the GSE relation, , making a rescaling , and after a lengthy calculation, we obtain the momentum distribution of 1D bosons up to order 1/c

where is the fermion two-particle correlation. Upto order 1/c, the normalization factor for strongly interacting 1D bosons has been calculated in Ref. [19] as .

By applying Wicks theorem and getting partial derivatives, the asymptotic behavior of the momentum distribution becomes

In the above calculation Eq. (42), we have proved that the leading order part gives the universal asymptotic behavior Eq. (43), whereas the second term in Eq. (42) is zero. According to the definition of the momentum distribution, equation (43) is equal to the relation given in Refs. [37] and [36]. Moreover, the Hellmannn–Feynman theorem relates to the derivative of the energy with respect to the interacting strength. So we have a two-particle local correlation function g₂(x₁,x₁ ) = n² e′(γ) where the derivative of total energy .^[66] The result of Eq. (43) shows that the tail of the momentum distribution decays with a power of 1/(p/ħ)⁴ and it has a sub-leading contribution determined by the GES parameter α.

A powerful universal relation connects the strength of short-range two-particle correlations to the thermodynamics of a many-particle system with zero-range interaction.^[55] The contact allows for a simple analytic expression that unveils how to scale the momentum distribution function at large momentum.^[54] We recall the Lemma discussed in Refs. [36], [46], and [68].

In the case of multiple singular points of the same type, the asymptotic behavior of the integral is given by the sum of all the corresponding contributions given by the right hand side of Eq. (45).

Considering discontinuity of the derivative wave function, equation (36) is given by

The structure factor is a property that defines how an ensemble of atoms scatters incident radiation. Experimentally, it is usually measured by two-photon Bragg scattering.^[42] Due to the translational invariance of the system, the density is constant,^[52] the static structure factor is defined as

In summary, we have presented a pedagogical study of various higher order local and non-local correlation functions of 1D Lieb–Liniger Bose gas. Using the Bethe wave function, we have studied the non-local 2M point correlation functions and the large momentum tail of the momentum distribution in term of the fractional statistical parameter α for the system with a strong repulsion. We have also discussed the higher-order correlation functions for the weakly interacting Bose gas via the Bethe ansatz wave function. It turns out that the leading order of the large-momentum tail is determined by the contact, which is related to the short-distance behavior of the two-body density matrix. The two-body density matrix has been given in term of the fractional statistical parameter α, giving the subleading order 1/c correction in the strong coupling regime, where the interaction strength c ≫1. Our results show that the fractional statistical parameter α attributes to the sub-leading contributions of the local, non-local correlation functions, two-body contact, static structure factor, etc. Our method is applicable to other exactly solvable models.